### Function used to examine presence of mistakes
### Last Revised and commented: 14 June 2016, JBS
### PLACE: Places a vector of members (y) in a space (x). Used in ENLARGE 
### ENLARGE: This is the actual game/model (requires MSM library for rtnorm function)

##########################################################################################

# Place function
# Input: x -- a vector of 0's equal to the length of the desired space
# 		 y -- a vector of the member locations in that space  
# 		 E.g. For a space of length ten x <- rep(0,10), and a membership vector y <- c(2,4,8), function outputs c(0,1,0,1,0,0,0,1,0,0)

place <- function(x,y) {
	x[as.numeric(names(table(y)))] <- table(y)
	return(x)
}

##########################################################################################
### Enlarge function

 
enlarge <- function(space,members,sdnmem=6,a=0.5,d=0.5,sdp=3,R) {
	
	output.list <- vector("list",length=R)
	
	# Use Place function to place initial members on the line
	space <- place(space,members)
	
	# Loop over each round
	for (r in 1:R) {
		output <- matrix(NA,nrow=length(members),ncol=10)

		# Determine location of pivotal member for decision-making in organization
		xpiv <- quantile(members, probs = d, type=2)
		
		# Draw the actual position of the new candidate for membership
		# Lower and upper truncate the distribution (best with space of 100) 
		# New members at the very extremes behave strangely
		# If located at the very extreme, it becomes impossible to draw perceptions on one side of the mean
		# rtnorm truncates here so members aren't right at boundary of the space.
		# Set to work with space of 100, would need to be altered if using different space
		nmem <- round(rtnorm(1,mean=mean(members),sd=sdnmem,lower=3,upper=98),digits=0)
		
		# Draw each existing members' perception of the candidate -- A vector equal to the number of existing members
		# Set to work with space of 100, would need to be altered if using different space
		nmemi <- round(rtnorm(length(members),mean=rep(nmem,length(members)),sd=abs(members-nmem)/sdp,lower=1,upper=100),0) 
		
		# replicate current member configuration for each existing member 
		hyp <- rep(list(space),length(members))
		hyp.1 <- rep(list(space),length(members))

		# Create holding bins for hypothetical new memberships for each existing member		
		hypmemlist <- vector("list", length(hyp))
		hypmemlist.1 <- vector("list", length(hyp.1))

		# For each existing member, create a hypothetical membership vector that reflects what the member perceives 
		# new membership to look like if candidate is admitted
	
		for (i in 1:length(hyp)) {
	
			hyp[[i]][nmemi[i]] <- hyp[[i]][nmemi[i]]+1
			hyp.1[[i]][nmem] <- hyp.1[[i]][nmem]+1
		
			hypmemlist[[i]] <- rep(which(hyp[[i]]>=1), hyp[[i]][which(hyp[[i]]>=1)])
			hypmemlist.1[[i]] <- rep(which(hyp.1[[i]]>=1), hyp.1[[i]][which(hyp.1[[i]]>=1)])

		}
	
		# Calculate the new hypothetical pivot for each existing member if the candidate were to be admitted.
    	xpiv.hyp <- unlist(lapply(hypmemlist,quantile, probs=d, type=2)) #
    	xpiv.hyp.1 <- unlist(lapply(hypmemlist.1,quantile, probs=d, type=2)) #
        
            
		# Create vector of votes in favor of admitting candidate 
		# Existing members vote in favor if the distance to current pivot is greater than the distance 
		# to the perceived hypothetical pivot if candidate is admitted.
		votes.perc <- as.numeric(abs(members-xpiv) > abs(members-xpiv.hyp))
		votes.real <- as.numeric(abs(members-xpiv) > abs(members-xpiv.hyp.1))
		
		output[,1] <- members
		output[,2] <- nmem
		output[,3] <- nmemi
		output[,4] <- xpiv
		output[,5] <- xpiv.hyp
		output[,6] <- xpiv.hyp.1
		output[,7] <- abs(members-nmem)
		output[,8] <- votes.perc
		output[,9] <- votes.real
		output[,10] <- as.numeric(votes.perc != votes.real)
		
		output.list[[r]] <- output
		
		}

# function returns the final membership vector after going through all of the rounds.
return(output.list)

}

library(msm)

sim1 <- enlarge(space = rep(0,100), members = c(44,44,44,44,46,47,48,53,54,55,56,56,56), R=1000)

outmat <- do.call(rbind,sim1)

mistakes <- outmat[,10]
dist <- outmat[,7]
dist2 <- dist^2

mod1 <- glm(mistakes~ dist+ dist2, family=binomial(link = "logit"))

mod2 <- lm(mistakes~ dist+ dist2)

pdf("mistakes.pdf")
plot(dist, jitter(mistakes), ylab="Mistake (Yes=1)", xlab="Distance", yaxt="n")
lines(sort(dist), 7*fitted(mod1)[order(dist)], lwd=2, col="red")
axis(2, at=c(0,1),labels=c(0,1))
dev.off()



